On numerical integration of coupled Korteweg-de Vries System 
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Abstract 

We introduce a numerical method for general coupled Korteweg-de Vries systems. The scheme is 
valid for solving Cauchy problems for arbitrary number of equations with arbitrary constant coeffi¬ 
cients. The numerical scheme takes its legality by proving its stability and convergence which gives 
the conditions and the appropriate choice of the grid sizes. The method is applied to Hirota-Satsuma 
(HS) system and compared with its known explicit solution investigating the influence of initial con¬ 
ditions and grid sizes on accuracy. We also illustrate the method to show the effects of constants with 
a transition to non-integrable cases. 


(Nl Introduction 

o 

-^Coupled Korteweg-de Vries system equations form a class of important nonlinear evolution systems. Its 
^importance comes (physically) from the wide application field it covers and (mathematically) from includ¬ 
ing both (weak) nonlinearity and third order derivatives (weak dispersion). It describes the interactions 
>-of long waves with different dispersion relations. Namely it is connected with most types of long waves 
^with weak dispersion {uj{k) —> 0, {k) 0), e.g. internal, acoustic and planetary waves in geophysical 

^hydro dynamics. 

C3 


It was introduced by A.Maxworthy, L.Redekopp and P.Weldman in studying the nonlinear atmo¬ 
sphere Rossby waves. R.Hirota and J.Satsuma 0 give single- and two-soliton solutions to some version of 
the system. R.Dodd and A.Fordy 0 found an L-A pair for Hirota-Satsuma equations. S.B.Leble derived 
the cKdV system for different hydrodynamical systems with explicit expressions for the nonlinear and dis¬ 
persion constants |^. He also developed the approach to the cKdV integration. S.B.Leble, S.P.Kshevetskii 
1^, ^ used the system in investigation of nonlinear internal gravity waves. A. Perelemova ^ used it in 
description of interaction of acoustic waves with opposite directions of propagation in liquids with bubbles. 


Others deal with integrability of the system (Dodd R. and Fordy A. [Q) from a Lax pair point of view, 
S.B.Leble 0 in Walquist-Estabrook theory. Foursov MV |@] described a new method for constructing 
integrable system of differential equations that reduced to cKdV equations. Oevel W. ^ considers inte- 
grable system of cKdV and found an inhnite hierarchy of commuting symmetries and conservation laws in 
involution. Zharkov A Yu. [^] obtained a new class of integrable KdV-like systems. Metin Curses, Ata- 
lay Karasu |]TI| found inhnitely many coupled system of KdV type equations which are integrable. They 
also give recursion operators. In studying the Painleve test classihcation of the system, Ayse (Kalkanli) 





Karasu [^| found new KdV systems that are completely integrable in the sense of WTC paper. He was 
looking for the integrable subclass of KdV systems given by Svinolupov |^. The later has introduced 
a class of integrable multicomponent KdV equations associated with Jordan algebras. John Weiss 
derived the associated ’’modihed” equations for HS system and from these the Lax pair is also derived. 
B.A.Kupershmidt |]^ showed that a dispersive system describing a vector multiplet interacting with the 
KdV held is a member of a bi-Hamiltonian integrable hierarchy. 


The signihcant achievement in numerical solution of the single KdV equation starts from the famous 
paper of Norman Zabusky and Martin Kruskal |^. It develops the idea of soliton solutions set for the 
integrable equations and enlighten the problems of effective integration scheme elaborating. The paper 
launched the numerous investigations and inventions in this held. Perhaps, the last publication that 
develop applications of recent theoretical achievements in numerical integration schemes is based on the 
notion of isospectral deformations [|^. Recently a multisymplectic twelve points scheme was produced 
T7| . This scheme is equivalent to the multisymplectic Preissmann scheme and is applied to solitary waves 


over long time interval. 

Shaohong Zhu |]^ had produced a diherence scheme for the periodic initial-boundary problem of the 
Coupled KdV (Hirota- Satsuma case) system. He use the inner product of the discrete function to obtain a 
scheme keeping two conserved quantities. His scheme is a nonlinear algebraic system for which a catch-ran 
iterative method is designed to solve it. 

The coupled KdV system representing most possible physical application ( related to the weak nonlinear 
dispersion ) to be considered in this work takes the following general form 


{dn)t + Cn{On)^ + '^gmkndkidm)^ + dn{9n)^^^ = 0 , U, m,/c = 1, 2, 3, ..., V, (1) 

k^m 

where 6n (x, t) is the amplitude of the wave mode as a function of space x and time t respectively. The 
constants c„ are the linear velocities and gmkn, are the nonlinear and dispersion coefficients. 


In the present work we introduce a numerical tool for solving coupled KdV system which is a develop¬ 
ment of the two step three time levels as Lax-Wendroff scheme p, ^ . Proving the theorem about stability 
and convergence of the scheme gives the opportunity to use it for different applications like Cauchy prob¬ 
lems for arbitrary number of equations and a wide class of initial conditions 9^ (x, 0). We consider in our 
problem an inhnite domain while the initial condition goes quickly enough to zero following the relation 

/(i + N)|«U.o)|* = o<oo 

keeping in mind the choice of smooth and integrable function. As an important corollary of the theorem 
one obtains conditions that have to be taken in account in choosing grid sizes. This numerical method is 
checked by applying it to HS system for which a good number of explicit solutions exist . We examine 
also the effects of equations coefficients and conditions of the problem on the solution. 

In the section 2 we introduce the difference scheme for arbitrary number of coupled KdV equations. 
We investigate stability and prove the convergence giving the condition have to be taken in account in 
choosing the grid sizes and how they are related. In section 3 we analyze HS system with two-parameters 
one-soliton explicit solution. The numerical method is applied to HS system and compared with the 
explicit solution. We analyzed the effects of the two parameters and initial condition on the form of the 
resulting solitons as will as on accuracy and show the results by hgures. We also produced (numerically) 
a multi-soliton solution for HS system and used the conservation law to estimate the expected number 
of solitons which agreed that we already obtained. Proving stability and convergence besides testing the 
results for HS system allows us in section 4 to use scheme for general cKdV system. Hence we illustrate 
























by plots the results of applying the scheme to slightly nonintegrable cKdV systems and other for a system 
with non-smooth initial conditions. 


2 The numerical method 

2.1 The difference scheme 

For the cKdV system (1) we introduce a numerical ( hnite-difference ) method of solution. A scheme 
which is two steps three time levels similar to the Lax-Wendroff one |^, The usual Lax-Wendroff is 
modihed such that the order of the hrst derivative becomes of order 0{Ax‘^). The approximation of the 
nonlinear terms is changed in such a manner that the integral of 6^ be a conserved one. The approach 
gives a solution that can be considered as some generalized solution, in the sense of Shwartz distribution 
theory, where the dispersion constants vanishes. This scheme is suitable to a nonlinear equations and is 
valid for n equations with arbitrary coefficients. The scheme can be simply derived beginning from Taylor 
series expansion as 

(9„)f ‘ = m + A« ((9,.),)J + 0|(Ai)=]. (2) 

where i and j are used to locate a point in the discrete domain and At is the time step while the subscript 
t means time derivative. Substituting for {On)^ in (2) using the system (1) to obtain 

mi*' = (mi - Ai |'c„ (9„), + 5 ; (9™), + d„ + O [(A i)"] . (3) 

The difference scheme is elaborated applying Lax idea for a half time step and leap frog method to the 
remaining half time step. In both steps {On)^ and are replaced by forth order O (Ax^) and second 

order accurate O (Ax^) central difference expressions. Hence (3) gives the following difference scheme 

{mA - (en)i) i\ +c„ /2ft +mi - m)Q /2ft 

k,m 

+ Cn ((6'n)-+2 ~ 2 (dn)i+i + 2 (dn)i_i “ (^n)i-2) = 0, e„ = (dn - , (4.1) 

where n, m, k are the modes numbers; i and j are discrete space and time respectively. The time step At is 
replaced for simplicity by t while h denotes spatial step. The equation (4.1) is accompanied with discrete 
equation for the intermediate layer as: 


mi*' - mi) n +c,. (mitt - mitt) / 2 ft+(o, 




mitt - (mitt) / 2 ft 


k.m 


+ e„ ( m)iti - 2 (mitt + 2 (mrnt - mu) tm = o (4.2) 




i+h 




2.2 Stability and convergence analysis 

For simplicity of the analysis we start by considering one equation of the system and give the details of 
stability and convergence. Then we apply the idea to the general cKdV system because it is rather close 
to that for one KdV equation but more bulky. 


2.2.1 Stability analysis for KdV scheme 

Consider one KdV equation of the system (1) 


+ 9 — 0 (5) 

Note again that the investigation we perform can be generalized for the case of any hnite number of modes. 
Considering the numerical scheme applied for the equation (5) 

- n)/r + - Ci)/2A + gOiW+i - Ci)/2A + e(»l+2 “ 2^1 + 2Ci “ «L2)/2A’ = 0 (6) 

Let us select a suitable norm. For this multiply equation equation (5) by 9 and integrate to yield 

d'^dx = 0 or / 6 ‘^dx = const, 

2 dt 

hence, by dehnition of L 2 norm,( 110112 )^ = J^^9‘^dx ,it may be written as (II^IU^ = const., i.e. the 
norm ||d ||2 is conserved and the equation can be treated in the L 2 norm. 

Now we will prove stability with respect to small perturbations ( because we consider nonlinear equa¬ 
tions ) of initial conditions. Strictly speaking it is the boundness of the discrete solution in terms of small 
perturbation of the initial data. So let us consider the differential 

. r = ..., 2 - 1 , 2 , 2 + 1 ,... (7) 

r 

f dOU 1 

for equation (5) denoting , dO^ = < dOj > and use WdO^l = ( {dOl)^ h 

d9U ^ 

dOl +2 > 

Rewrite also the relation (7) in the matrix form 

r 

where dO^ is a small perturbation of initial data and the subscripts are omitted for simplicity. 

Stability requires the boundedness of the product T*" in a sense that the norm is bounded 

r r 

by some constant, i.e. < C. Here C is a constant, and the matrix norm is a spectral norm. 

r 

For this the sufficient condition is ||T’^|| < where a is a constant, that is independent of r. The 
case ||T’’|| < jg ^ sufficient condition of stability also, but only if |a(r, h)| < const < 00 . If 

|a(r,/i)| < const, including r, h ^ 0, for some dependence r = f{h), then we can say about conditional 
stability. Namely this kind of stability will be climbed below. 

To calculate T^, rewrite the scheme ( 6 ) in the form ( 6 *^+ 2 , ^i+i, dj, di_^, 6 i_ 2 ). 

So 

= 5,,, - (cr/2h) - {gT/2h) \9i + 5,,, {91, - 9l_,)] 

— {cT/2h?) [5ij^2,r ~ 25j+l,r + — 5i-2,r\ ■ (8) 




Rewriting (8) in terms of the identity (E), symmetric (S) and anti-symmetric (A) matrices 
T^'+i = C + 


i(0i - «l+i) ■S.+ix - (n - C.) i-i, + [9|+1 - 9'_i]) 


{Aj+y., = ~ ((«;+i.+., - (ei + 




2 h 

er 


Ah 

[A-|-2,r -|- A—2,r] 


ll^^+^ll < \g\T max(l0^ .[ , [0'^^ .j) , /(h) and • = [61^^ - /(2/i) 


one arrive at 




h 


h 




||T^'+i||^ = II (T^+i)*T^+i|| = II (E - + S^+^) {E + + S^+^) || 

< l + 2||^^’+i|| + (||A^’+i|| + ll^^'+^ID" 

< 1 + 2 1^1 r max |dR| + ^'2 \g\ max \eE\ + max | + y ^ 


< 6“"^ where a = 2 |5(| max Id^^-I-h r ( 2 |5(| max Id^^jl-Fmax l^^ l-h 


which is a necessary condition of stability. The scheme is stable if a < constant < cx) in spite of r, h —0. 
This is a conditional stability of the scheme. It means that it is required for stability that r —*> 0 more 
faster than h ^ 0, or 


r < (constant), h®, constant < 00 (9) 

Therefore, for small enough r we can simplify the expression for a 

a = 2\g\ max|d^_.|+r 

In practical calculations the time step r should be chosen so that it would satisfy r = 0(1), where 

to is the time of simulation (0 < t < to)- In future, when we shall be suggesting some better numerical 
scheme, we will essentially use our observation that stability depends only on the dispersion terms. And 
now we will try to accomplish our short investigation of the scheme (6) by strong proving of the numerical 
scheme convergence. 

2.2.2 The proof of the KdV scheme convergence 

Now we prove that a solution of equation (6) converges to a solution of (5), if the exact solution is a 
continuously-differentiable one. Let us denote by 9(x,t) a solution of the equation (5). We substitute 6*^ = 
9(xi,tj) + vj into (6), is a error between the difference solution 9^ and the exact solution 9(xi,tj)- We 




obtain the equation for 

+ c {vy - vy) / 2 h + ge{xi, tj) {vy - vy) / 2 h + gvj {e{xi+i,tj) - e{xi_i,tj)) / 2 h 
+ gvj {vy - vy) / 2 h + e {vy - + M -1 - vi- 2 ) /2h^ = 

- ((9(xi, tj+i) - e{xi, tj)) /t + c {e{xi+i,tj) - 9{xi_i, tj)) /2h + g9{xi, tj) {9{xi+i,tj) - 9{xi_i,tj)) /2/i 

+e ( 9 (xi+ 2 , tj) - 29{xi+i,tj) + 29{xi_i, tj) - 9 {xi_ 2 , tj)) / 2 h^) 

Let us take into account that 

vj - T {c{vy - v\y/ 2 h + g9{xi, tj){vy - v\y/ 2 h + gv\{9{xi+^,tj) - 9{xi_x, tj))/ 2 h 

+e{vy - 2vy + 2vy - vy)/2h^) = ^ 

k 

Using the operator T-^+i introduced above, this equation may be rewritten in the form 




I + 9vl {vy - viy I2h = 


- ((9(xi, tj+i) - 9(xi, tj)) /t + c {9{xi+i,tj) - 9{xi_i,tj)) /2h + g9{xi, tj) {9{xi+i,tj) - 9{xi_i,tj)) /2h 

+e { 9 {xi+ 2 , tj) - 29{xi+i,tj) + 29{xi_i,tj) - 9 {xi_ 2 , tj)) / 2 h^) 
The right part of this relation is a quantity of order 0(r + h?). So, we can write 




jr + gvl {vy - viy / 2 h = 0 {t + h?), 


or 


^ _ rfi, /f = - 0(r + ft^). 


v, = 


2 h 


We hnally arrive at the inequality that compare the norms: 

■^11 < 11'^'^ iT 

/l2 

To explain how this estimate was obtained, follow the expressions 


( 10 ) 


U - 


Yinfh] 



h +0(r + /i^) 


^ 1^1 h 




+ 0{T + h^). (11) 


Using the Schwartz inequality \\AB\\ < ||A|| ||-B||, the formulas (10) could be transformed as 

l|,„i+i|l < ||r,+i|| ||-(,^|| + T||y|| < ||T^+‘|| iiT-'ii + T(||r+‘|| ||f-‘|i + iifii) < 

\\T>*'\\ ||T^|| ||P-‘|| +T(||r+‘[| [|r|| llfll + IjP+'ll Ilf-‘II + Ilf II < 

6“'^ |k.“|| +r(e“'<‘-‘' ||/“|| ||/‘|| + ...||f ||) < 


||t;°l| +Mmax(||/^||) < 

k<j 


oO-rj 1L,0| 


+ M (^M ||«‘+‘ f + 0(t + ftf j , M = 


( 12 ) 



To derive (12), we have used the iteration of the hrst of formula (we substituted the formula into itself, 


but for index less than 1) and using (11). Then we have utilized the formula for a sum of geometric series. 
Farther the inequality obtained in (12) has a solution 


„i+i| 


< I 1 - W 1 - 4m|j (6“-^ ||u0|| + MO{t + h2)) I / 


(2mM 

V h'2 


If we take into account (9), and use||n°|| = 0,we obtain 


-2+MI < 1- 


|»W+/,=))/(2^ |»| 


hi 


hi 




hi 


'-MO{t + h") / 


)'e 


hi 


= M * 0{t + h^) 


The constant M is bounded, in spite of j —> cx), because of jr < oo. Therefore, the convergence is proved. 

2.2.3 The coupled KdV scheme 

The numerical scheme for the system of the cKdV (4.1)-(4.2) is also conditionally stable and convergent 
one. The proof for this scheme is close to that given before, but a bit bulky. We deal with a vector: 

t/={hi 92 ... hiv y 

as a dependent variable instead of the simple variable 9 in the case of one KdV. For this vector case the 
norm used has the form 

Ill'll = 

\ i=l i 

The conditions connecting time and space steps for this scheme look also similar but with different constants 

81 * 

max(|e„|) * * to = 0(1). 

n 4* h^^ 


h 


3 Checking the numerical method 

The numerical method is tested by applying it to Hirota-Satsuma system. Namely the two parameters 
one-soliton explicit solution is used . 


3.1 Analytic solution ( explicit formula ) of Hirota-Satsuma system 

Darboux transformation (DT) that account a deep reduction for this specihc HS case of cKdV is used 


^ to obtain explicit solutions to HS system. The Lax representation of the HS equations is based on 
the matrix 2x2 spectral problem of the second order. For this problem the deep reduction scheme is 


applied (with the help of the conserved bilinear - forms) and supports the constrains on the potential while 
the iterated DT are performed. The iterated DT in determinant form and the covariance of the bilinear 
forms with respect to DT under restrictions gives N soliton solution of HS system. The system of HS we 
use to check the scheme has the form: 


(»,), - 0.26 (»,)3. - 1.5 {0,1 (0,) + 3 (02)„ (02) = 0 


(13.a) 


















(13.b) 


(02\ +0.5 (^ 2 ) 3 . +1.5 {92), (0i) =0. 

This system has two-parameters one soliton solution: 

61 = —2m^(—1 + d'^ + 2d Sin{Xi) * Sinh{\ 2 ))/{d Cos(Ai) + Cosh{X 2 ))'^ , 

02 = {2 + 2d^y^my{dCos{Xi) + Cosh{X2)), (14) 

Ai = .5mH + mx and X 2 = .5m^t — mx. 

with real constants m, d. For small \d\ this solution is a smooth function but for \d\ > 1 poles appear .The 
following hgures show some choice of m and d to show the effect of these two parameters on the solution. 
Figure 1 show that, for constant d, the amplitude is proportional to m while the wave width is inversely 
proportional to it. Figure 2 show that, for the given m, d affects on soliton shape, namely the hrst mode, 
while the amplitude of the second is inversely proportional to d. 




Fig.l. For a constant d, the amplitude is proportional to m while the wave width is inversely 

proportional to m. 




Fig. 2. For the same m, d affects on soliton shape, namely the hrst mode, while the amplitude 

of the second mode is inversely proportional to d. 


3.2 Calculations by numerical scheme and comparison results 

HS system (13) is solved numerically using the scheme (4) with initial condition from formula (14) (t = 0) 
and the results are compared with the explicit solution. It is found that, keeping the restriction on the 
choice of r and h and relation between them, the initial wave modes amplitude affects on the accuracy 
of the results . Also the error decreases as the mesh is rehned. Namely smaller amplitude ( of order one ) 
gives better results as shown below in hgure (3). It gives the percentage error calculated as follow 


%Error 


\Explicit solution — Numerical solution\ 
Initial amplitude 


We relate the error to the initial amplitude to show the physical signihcance of the error. 











Fig.3 % error is proportional to the amplitude (A). 

The plots show that the error is proportional to the amplitude (A), where as shown in figure the maximum 
relative error in the case (A= 2 ) is 1 % while in the case (A=3.4) is 3 %. The reason may be due to the 
higher velocity in the larger one, hence more interactions impact. It also shows that the error increases near 
the peak points. The reason of these oscillations in plots appearance is that the numerical and analytical 
plots intersect over the space domain. 

4 Applying the scheme to different applications 

Stability analysis and checking performed to the scheme in the general cKdV equations give the ability 
of using this scheme to solve other problems for which analytical solutions have not been found. We first 
consider the multi-soliton solution decay of the localized initial condition for the single KdV equation of 
HS system (figure 4.a) then for the complete HS system (figure 4.b). In both we use the initial condition 
from formula (14) but with 10 times the width and twice the amplitude. 



Fig. 4.A The multi-soliton decaying of the isolated (first) KdV of HS system (13.a) 



Fig. 4.b The multi-soliton decaying of the complete HS system (13). 

The second mode affects (interacts) the first one which results in the right direction soliton-like ’’tail” as 
shown in the figure 4.b. We estimate, using the conservation low (derived below), the expected number of 
solitons that already obtained in the numerical solution as 

01 X (13.1) -202 X (13.2) ^ 

-^[■59f — 62 ] + — • 2501013 , 3 , + • 1250 ^ 3 , — 9202xx + - 5013 ,] = 0 

01 and 02 —> 0 as a: —> ±cx) 

/oo°° (-^^1 ~ 9l)dx = const. 
















































Next we go to the solution of nonintegrable HS system. The integrable HS system (13.a,b) may be shifted 
to ’’slightly” nonintegrable one by small change of the dispersion constant of the hrst equation to have the 
new nonintegrable HS system 

(0i), - 0.2 (0i)3, - 1.5 (0i), (0i) + 3 ( 62 ), ( 62 ) = 0 (15.a) 

(^ 2 )^ + 0.5 ( 02 )^^ + 1.5 ( 62 )^ (Oi) = 0. (15.b) 

Using our scheme with initial condition from (14) we hnd that the scheme works satisfactory (in the sence 
of convergence) even for nonintegrable HS system as shown from Figure 5 below. The solution looks like 
a soliton one for small time. 



Fig. (5) The numerical solution of integrable and slightly nonintegrable HS system. 

Also the solution using a non-smooth initial condition for HS system (13) is shown in hgure 6 below . 



Fig.6 The numerical solution of Non-smooth initial condition for HS system (13). 
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